##### read files ##### -----
rm(list=ls())
setwd("D:/OneDrive - zju.edu.cn/data/public_database/cBioportal/brca_tcga")

metadata_patient = read.csv("data_bcr_clinical_data_patient.txt",skip=4,header = 5,sep="\t")
metadata_sample = read.csv("data_bcr_clinical_data_sample.txt",skip=4,header = 5,sep="\t")
exprset_median = read.csv("data_RNA_Seq_v2_expression_median.txt",header = T,sep="\t")
exprset_zscore = read.csv("data_RNA_Seq_v2_mRNA_median_Zscores.txt",header=T,sep="\t")

##### data format ##### -----
choose_columns_patient=c("PATIENT_ID","ER_STATUS_BY_IHC","PR_STATUS_BY_IHC","IHC_HER2","OS_STATUS","OS_MONTHS","DFS_STATUS","DFS_MONTHS")
choose_columns_sample=c("PATIENT_ID","SAMPLE_ID")
metadata_patient = metadata_patient[,choose_columns_patient]
tnbc_patient = metadata_patient[metadata_patient$ER_STATUS_BY_IHC=="Negative"&metadata_patient$PR_STATUS_BY_IHC=="Negative"&metadata_patient$IHC_HER2 =="Negative",]
metadata_sample = metadata_sample[,choose_columns_sample]
rownames(tnbc_patient) = tnbc_patient[,1]
tnbc_sample = merge(tnbc_patient,metadata_sample,by="PATIENT_ID")
rownames(tnbc_sample) = as.character(tnbc_sample$SAMPLE_ID)
rownames(tnbc_sample) = gsub("\\-", ".", rownames(tnbc_sample), ignore.case = FALSE, perl = FALSE,fixed = FALSE, useBytes = FALSE)
tnbc_sample = tnbc_sample[-48,]
write.csv(tnbc_sample,"tnbc_sample.csv")
tnbc_sample = read.csv("tnbc_sample.csv",header = T,row.names = 1)

matrix_median = as.matrix(exprset_median[,3:ncol(exprset_median)])
rownames(matrix_median) = exprset_median[,1]
tmp <- by(matrix_median,
          rownames(matrix_median),
          function(x)x[which.max(rowMeans(x)),])
matrix_median <- do.call(rbind,tmp)

tnbc_median = matrix_median[,rownames(tnbc_sample)]

###
matrix_zscore = as.matrix(exprset_zscore[,3:ncol(exprset_zscore)])
rownames(matrix_zscore) = exprset_zscore[,1]
tmp <- by(matrix_zscore,
          rownames(matrix_zscore),
          function(x)x[which.max(rowMeans(x)),])
matrix_zscore <- do.call(rbind,tmp)

tnbc_zscore = matrix_zscore[,rownames(tnbc_sample)]

write.csv(matrix_median,"matrix_median.csv")
write.csv(matrix_zscore,"matrix_zscore.csv")
write.csv(tnbc_median,"tnbc_median.csv")
write.csv(tnbc_zscore,"tnbc_zscore.csv")

##### survival ##### -----
rm(list=ls())

library(survival)
library(survminer)
library(GSVA)

tnbc_median = read.csv("tnbc_median.csv",header=T,row.names = 1)
tnbc_zscore = read.csv("tnbc_zscore.csv",header=T,row.names = 1)
tnbc_sample = read.csv("tnbc_sample.csv",header = T,row.names = 1)
tnbc_sample = tnbc_sample[,5:8]

matrix = as.matrix(tnbc_median)

IFNA_gene = data.frame(gene=c("HLA-G","PTPN2","IRF9","LGALS3BP","NUP210","EIF2AK2","KPNA7","CXCL10","NLRC5","UBE2L6",
                              "B2M","IFI6","IL-15","IRF1","STAT1","C1S",
                              "GBP6","GMPR","USP18","LAMP3"))
IFNA_score = gsva(matrix,IFNA_gene)

high_score = as.data.frame(IFNA_score[,IFNA_score["gene",]>0])
colnames(high_score) = "score"
IFNA_high = tnbc_sample[rownames(high_score),]
IFNA_high$group = "high"
IFNA_high = cbind(IFNA_high,high_score)
IFNA_high = IFNA_high[order(IFNA_high$score),]

low_score = as.data.frame(IFNA_score[,IFNA_score["gene",]<0])
colnames(low_score) = "score"
IFNA_low = tnbc_sample[rownames(low_score),]
IFNA_low$group = "low"
IFNA_low = cbind(IFNA_low,low_score)
IFNA_low = IFNA_low[order(IFNA_low$score),]

dat = rbind(IFNA_high,IFNA_low)

my.surv <- Surv(dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')
kmfit1 <- survfit(my.surv~group,data = dat)
curve = ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE,legend.title="",xlab = "Months") +
  ggtitle("IFNA score survival curve")

pdf("./tcga IFNA score.pdf", onefile = FALSE)
print(curve)
dev.off()
curve



